Assignment 02

Name :- ADIKARI A.M.A.D

Index No :- 190021A

# import libraries
import numpy as np
from scipy.optimize import minimize
from scipy import linalg
import matplotlib . pyplot as plt
import cv2 as cv

Question 01

# function for detecting inlier points
def inliers(points, center,  r, thresh_dis):
    low_thresh, high_thresh = r-thresh_dis, r+thresh_dis   # threshold range
    dis = np.sqrt((points[:, 0]-center[0])**2+(points[:, 1]-center[1])**2)  # calculate the distance from center to points
    indices = [i for i in range(0, len(points)) if low_thresh<=abs(dis[i])<=high_thresh]        # indices of inlier points
    inl = np.array([list(points[i]) for i in indices])      # inlier points
    return inl

# function for getting points of circles
def circle_points(center, r):
    theta = np.linspace(0, 2*np.pi, 100)
    # cordinate points of the circle
    x = center[0] + r*np.cos(theta)
    y = center[1] + r*np.sin(theta)
    return (x, y)

# function for calculating the mean absolute error
def mean_absolute_error(points, center, r):
    g, f = center[0], center[1]
    dis = np.sqrt((points[:, 0]-g)**2+(points[:, 1]-f)**2)
    error = np.sum(abs(dis-r))/len(points)
    return error

# function for finding best RANSAC fitting with maximum inlier count
def best_RANSAC(ran):
    max_inl_count = 0
    best_ran = [[[], []], [[], []], [], [], [[], []]]
    for i in range(0, len(ran[0][0])):
        if max_inl_count<len(ran[2][i]):
            max_inl_count = len(ran[2][i])  # maximum inlier count
            for j in range(0, 5):
                if j==2 or j==3:
                    best_ran[j] = ran[j][i]
                else:
                    best_ran[j][0] = ran[j][0][i]
                    best_ran[j][1] = ran[j][1][i]
        # if inlier counts are equal then consider the minimum mean absolute error
        elif max_inl_count==len(ran[2][i]):
            if (best_ran[3] > ran[3][i]):
                max_inl_count = len(ran[2][i])
            for j in range(0, 5):
                if j==2 or j==3:
                    best_ran[j] = ran[j][i]
                else:
                    best_ran[j][0] = ran[j][0][i]
                    best_ran[j][1] = ran[j][1][i]

    print('inliers count :',max_inl_count)
    return best_ran

# function for getting RANSAC circles
def RANSAC(points, iterations, inlier_count, max_thresh, cand=0):
    max_r = (abs(points[:, 0].min())+abs(points[:, 0].max()))/2+5  # max radius
    ran = [[[], []], [[], []], [], [], [[], []]]    #format : center-RANSAC, center-best, sample-RANSAC, sample-RANSAC, inliers-RANSAC, error-RANSAC, radius-RANSAC, radius-best
    for i in range(0, iterations):
        p1, p2, p3 = points[np.random.choice(points.shape[0], 3, replace=False), :]  # 3 random points
        
        # calculate the center and the radius using matrices
        A = np.array([[2*p1[0], 2*p1[1], 1], [2*p2[0], 2*p2[1], 1], [2*p3[0], 2*p3[1], 1]])
        B = np.array([[p1[0]**2+p1[1]**2], [p2[0]**2+p2[1]**2], [p3[0]**2+p3[1]**2]])
        res = np.linalg.pinv(A) @ B
        g, f = res[0][0], res[1][0]
        r = np.sqrt(res[2]+g**2+f**2)
        # valid radius (below the max radius)
        if (r[0]>max_r):
            continue
        inlier_points = inliers(points, (g, f), r, max_thresh) # inlier points

        # check that inlier counts are greater than or equal to the threshold value
        if (inlier_count<=len(inlier_points)):
            # cand=1 : new candidate circle for a higher accuracy
            if (cand==1):
                # RANSAC circle datas
                ran[0][0].append((g,f))
                ran[1][0].append(np.array([[p1[0], p1[1]], [p2[0], p2[1]], [p3[0], p3[1]]]))
                ran[2].append(inlier_points)
                error = mean_absolute_error(inlier_points, (g, f), r)
                ran[3].append(error)
                ran[4][0].append(r)
            # cand = 0 : first candidate circle for the selected random sample
            if (cand==0):
                # datas for the best sample
                ran[0][1].append((g, f))
                ran[1][1].append(np.array([[p1[0], p1[1]], [p2[0], p2[1]], [p3[0], p3[1]]]))               
                ran[4][1].append(r)
                ran2 = RANSAC(inlier_points, iterations, inlier_count, max_thresh, 1)
                for i in range(0, len(ran2[0][0])):
                    ran[0][0].append(ran2[0][0][i])
                    ran[0][1].append((g, f))
                    ran[1][0].append(ran2[1][0][i])
                    ran[1][1].append(np.array([[p1[0], p1[1]], [p2[0], p2[1]], [p3[0], p3[1]]]))
                    ran[2].append(ran2[2][i])
                    ran[3].append(ran2[3][i])
                    ran[4][0].append(ran2[4][0][i])
                    ran[4][1].append(r)
    return ran

# sample points
# np . random . seed ( 0 )
N = 100
half_n = N// 2
r = 10
s = r/16
t = np.random.uniform(0, 2*np.pi, half_n)
n = s*np.random.randn(half_n)
x, y = (r + n)*np.cos(t), (r + n)*np.sin(t)
X_circ = np.hstack((x.reshape(half_n, 1), y.reshape(half_n, 1)))
m, b = -1, 2
x = np.linspace(-12, 12, half_n)
y = m*x + b + s*np.random.randn(half_n)
X_line = np.hstack((x.reshape(half_n, 1), y.reshape(half_n, 1)))
X = np.vstack((X_circ, X_line))

X = X[np.argsort(X[:, 0])] # sort


s = 3       # points for a sample
e = 0.6     # outlier ratio
p = 0.99    # probability of success
ite = int (np.log(1-p)/np.log(1-(1-e)**s)) # no. of iterations
thresh = 1  # distance threshold
inliers_count = N*(1-e) # inliers count for a circle
print('number of iterations :', ite)

# getting RANSAC and best circle fitting
data = RANSAC(X, ite, inliers_count, thresh)
best_ran = best_RANSAC(data)

center_ran, center_best = best_ran[0][0], best_ran[0][1]
sample_ran, sample_best = best_ran[1][0], best_ran[1][1]
inlier_points = best_ran[2]
error = best_ran[3]
r_ran, r_best = best_ran[4][0], best_ran[4][1]
circle_ran = circle_points(center_ran, r_ran)
circle_best = circle_points(center_best, r_best)

fig, ax = plt.subplots(figsize=(15, 15))
ax.plot(X[:, 0], X[:, 1], 'bo', label='points')
ax.plot(center_ran[0], center_ran[1], 'm+', markersize=10)
ax.plot(circle_ran[0], circle_ran[1], 'm', label='RANSAC')
ax.plot(center_best[0], center_best[1], 'b+', markersize=10)
ax.plot(circle_best[0], circle_best[1], 'b', label='Best sample')
ax.plot(inlier_points[:, 0], inlier_points[:, 1], 'yo', label='Inliers')
ax.plot(sample_best[:, 0], sample_best[:, 1], 'ro', label='Best sample')
ax.set_aspect('equal')
ax.legend(loc='best', fontsize=20)
ax.set_xlabel('x', fontsize=15)
ax.set_ylabel('y', fontsize=15)

plt.show()

number of iterations : 69
inliers count : 44

Question 02

# function for getting coordiantes of clicked points
def mouse_click(event, x, y, flags, params):
    global counter
    # check if the left button of the mouse is clicked
    if event == cv.EVENT_LBUTTONDOWN:
        clicked_points[counter] = (x,y)
        counter += 1

def mousePoints(img):
    while True:
        for i in range(0, 4):
            # show a circled on clicked points
            cv.circle(img, (clicked_points[i][0], clicked_points[i][1]), 2, (0, 0, 255), cv.FILLED)
        
        if counter==4:
            cv.destroyAllWindows()
            break
        
        # showing the image
        cv.imshow('image', img)
        cv.setMouseCallback('image', mouse_click)
        cv.waitKey(1)

def warpAndBlend(img1, img2, alpha, beta, gamma):
    global counter
    global clicked_points
    clicked_points = np.array([(0, 0), (0, 0), (0, 0), (0, 0)]) # format : left_top, right_top, right_bottom, left_bottom
    counter=0
    img1 = cv.resize(img1, (1500, 900))
    # click points on the img1 image
    counter = 0

    img1_copy = img1.copy()
    mousePoints(img1_copy)

    # color converting
    img1 = cv.cvtColor(img1, cv.COLOR_BGR2RGB)
    img2 = cv.cvtColor(img2, cv.COLOR_BGR2RGB)

    img2_points = np.array([(0, 0), (img2.shape[1], 0), (img2.shape[1], img2.shape[0]), (0, img2.shape[0])])   # corner points of the img2

    h, status = cv.findHomography(img2_points, clicked_points)              # homography
    output_img = cv.warpPerspective(img2, h, (img1.shape[1], img1.shape[0]))        # warping img2 into img1
    blend_img = cv.addWeighted(img1, alpha, output_img, beta, gamma)           # blending img2 into img1

    fig, ax = plt.subplots(1, 2, figsize=(15, 8))
    ax[0].imshow(img1)
    ax[0].set_xticks([])
    ax[0].set_yticks([])
    ax[0].set_title('main image', color='white')

    ax[1].imshow(img2)
    #ax[0][1].plot(img2_points[:, 0], img2_points[:, 1], 'go')
    ax[1].set_xticks([])
    ax[1].set_yticks([])
    ax[1].set_title('warping image', color='white')

    fig, ax = plt.subplots(figsize=(10, 10))
    ax.imshow(output_img)
    ax.set_xticks([])
    ax.set_yticks([])
    ax.set_title('Warped image', color='white')

    fig, ax = plt.subplots(figsize=(25, 15))
    ax.imshow(blend_img)
    ax.set_xticks([])
    ax.set_yticks([])
    ax.set_title('Blended image', color='white')

    plt.show()

# building and the flag
building = cv.imread('./Q2_images/001.jpg')
flag = cv.imread('./Q2_images/flag.png')
assert building is not None
assert flag is not None
warpAndBlend(building, flag, 1, 0.8, 0)

# car and the logo
car = cv.imread('./Q2_images/ferrari.jpg')
logo = cv.imread('./Q2_images/horse_logo.jpg')
assert car is not None
assert logo is not None
# adding a mask
color = [200, 0, 0]
mask = cv.inRange(logo, np.array([0, 0, 0]), np.array([150, 150, 150]))
logo[mask>0] = color[::-1]
logo[mask==0] = 0
warpAndBlend(car, logo, 0.7, 0.3, 0)

# wall and the headset
wall = cv.imread('./Q2_images/studio.jpg')
headset = cv.imread('./Q2_images/headset.jpg')
assert wall is not None
assert headset is not None
# adding a mask
mask = cv.inRange(headset, np.array([0, 0, 0]), np.array([127, 127, 127]))
headset[mask>0] = 255
headset[mask==0] = 0
warpAndBlend(wall, headset, 1, 1, 0)

Question 03

# function for getting matching pairs using SIFT algorithm
def sift_matching(img1, img2, im1, im2):
    # gray images
    img1_gray = cv.cvtColor(img1, cv.COLOR_BGR2GRAY)
    img2_gray = cv.cvtColor(img2, cv.COLOR_BGR2GRAY)

    # SIFT
    sift = cv.SIFT_create()
    bf = cv.BFMatcher(cv.NORM_L2SQR, crossCheck=True)
    key1, des1 = sift.detectAndCompute(img1_gray, None)
    key2, des2 = sift.detectAndCompute(img2_gray, None)
    matches = bf.match(des1, des2)
    matches = sorted(matches, key=lambda x:x.distance)
    # output image
    out_img = cv.drawMatches(img1, key1, img2, key2, matches[:100], img2, flags=2)

    print('matching pairs for img{}to{}'.format(im1, im2))
    fig, ax = plt.subplots(figsize=(20, 20))
    ax.imshow(out_img)
    ax.set_xticks([])
    ax.set_yticks([])
    plt.show()
    
    return matches, key1, key2

# function for abstracting the points of query and train image from the matching pairs
def matched_coordinates(matches, key1, key2):
    p1, p2 = [], []
    for match in matches:
        k1_index = match.queryIdx
        k2_index = match.trainIdx
        point1 = key1[k1_index].pt
        point2 = key2[k2_index].pt
        p1.append(point1)
        p2.append(point2)
    return np.array(p1, dtype=np.float32), np.array(p2, dtype=np.float32)

# function for count the inliers
def inliers(H, points1, points2, thresh_dis):
    count = 0
    inl = [] # inlier points
    for i in range(0, len(points1)):  
        # apply homography for every query image point and check the validity     
        B = np.array([[points2[i][0]], [points2[i][1]], [1]])
        A = np.array(H) @ B         
        Ax, Ay = A[0][0]/A[-1][0], A[1][0]/A[-1][0]
        Px, Py = points1[i][0], points1[i][1]
        if (np.sqrt((Ax-Px)**2 + (Ay-Py)**2)<=thresh_dis):
            count += 1          
            inl.append(i)
    return [count-4, inl]
    

# function for RANSAC algorithm
def RANSAC(matches, key1, key2,  iterations, thresh_dis):
    best_points = [[], []]
    max_inliers = 0
    indices = []
    best_H = []
    points1, points2 = matched_coordinates(matches, key1, key2)
    for i in range(0, iterations):
        # random points
        k1, k2, k3, k4 = np.random.randint(0, len(matches)-1, 4)
        q_index1, t_index1 = matches[k1].queryIdx, matches[k1].trainIdx
        point_q1, point_t1 = key1[q_index1].pt, key2[t_index1].pt
        q_index2, t_index2 = matches[k2].queryIdx, matches[k2].trainIdx
        point_q2, point_t2 = key1[q_index2].pt, key2[t_index2].pt
        q_index3, t_index3 = matches[k3].queryIdx, matches[k3].trainIdx
        point_q3, point_t3 = key1[q_index3].pt, key2[t_index3].pt
        q_index4, t_index4 = matches[k4].queryIdx, matches[k4].trainIdx
        point_q4, point_t4 = key1[q_index4].pt, key2[t_index4].pt

        q_points = [point_q1, point_q2, point_q3, point_q4] # query image points
        t_points = [point_t1, point_t2, point_t3, point_t4] # train image points

        # check the colinearity
        colinear =  bool(np.cross(point_q1, point_q2)) and bool(np.cross(point_q2, point_q3)) and bool(np.cross(point_q3, point_q4))
        colinear =  colinear and bool(np.cross(point_q2, point_q3)) and bool(np.cross(point_q2, point_q4))
        colinear =  colinear and bool(np.cross(point_q3, point_q4))
        if (not colinear):
            continue
        
        # form the A and B matrices
        A = []
        for i in range(0, 5):
            if i!=4:
                A.append([-t_points[i][0], -t_points[i][1], -1, 0, 0, 0, t_points[i][0]*q_points[i][0], t_points[i][1]*q_points[i][0], q_points[i][0]])
                A.append([0, 0, 0, -t_points[i][0], -t_points[i][1], -1, t_points[i][0]*q_points[i][1], t_points[i][1]*q_points[i][1], q_points[i][1]])
            else:
                A.append([0, 0, 0, 0, 0, 0, 0, 0, 1])
        A = np.array(A)
        B = np.array([[0], [0], [0], [0], [0], [0], [0], [0], [1]])

        # calculating the homography using A and B matrices
        h = np.linalg.inv(A.T @ A) @ A.T @ B
        H = [[], [], []]
        for i in range(0, 9):
            H[i//3].append(h[i][0])

        inliers_n, ind = inliers(H, points1, points2, thresh_dis) # get the inliers

        # select the sample with the maximum number of inliers
        if inliers_n>max_inliers:        
            best_points[0] = [point_q1, point_q2, point_q3, point_q4]
            best_points[1] = [point_t1, point_t2, point_t3, point_t4]
            max_inliers=inliers_n
            best_H = H
            indices = ind
            
    return best_H, [np.array(best_points[0], dtype=np.float32), np.array(best_points[1], dtype=np.float32)], indices

# function for implementing the graphs with data points and returning the homography
def homoAndGraph(matches, img1, img2, key1, key2):
    points1, points2 = matched_coordinates(matches, key1, key2) # query and train image points

    s = 4       # points for a sample
    e = 0.7     # outlier ratio
    p = 0.99    # probability of success
    ite = int (np.log(1-p)/np.log(1-(1-e)**s)) # no. of iterations
    thresh = 5  # distance threshold

    H, p, indices = RANSAC(matches, key1, key2, ite, thresh)
    inliers1, inliers2 = points1[indices], points2[indices]
    img1_h_points, img2_h_points = p[0], p[1] # sample points used to get the best homography
    
    print('matches count :', len(matches))
    print("iteration :", ite)
    print('inliers count :', len(inliers1))

    fig, ax = plt.subplots(1, 2, figsize=(15, 8))
    # data points spread of the query image
    ax[0].plot(points1[:, 0], points1[:, 1], 'bo')  # outliers
    ax[0].plot(inliers1[:, 0], inliers1[:, 1], 'yo')    # inliers
    ax[0].plot(img1_h_points[:, 0], img1_h_points[:, 1], 'ro')  # best sample points

    # data points spread of the train image
    ax[1].plot(points2[:, 0], points2[:, 1], 'bo')  # outliers
    ax[1].plot(inliers2[:, 0], inliers2[:, 1], 'yo')    # inliers
    ax[1].plot(img2_h_points[:, 0], img2_h_points[:, 1], 'ro') # best sample points

    fig, ax = plt.subplots(1, 2, figsize=(15, 8))
    # query image
    ax[0].imshow(img1)
    ax[0].plot(img1_h_points[:, 0], img1_h_points[:, 1], 'ro', markersize=5)    # best sample points
    # train image
    ax[1].imshow(img2)
    ax[1].plot(img2_h_points[:, 0], img2_h_points[:, 1], 'ro', markersize=5)    # best sample points

    plt.show()
    
    return H

# reading the images
img1 = cv.imread('./Q3_images/img1.ppm')
img2 = cv.imread('./Q3_images/img2.ppm')
img3 = cv.imread('./Q3_images/img3.ppm')
img4 = cv.imread('./Q3_images/img4.ppm')
img5 = cv.imread('./Q3_images/img5.ppm')
assert img1 is not None
assert img2 is not None
assert img3 is not None
assert img4 is not None
assert img5 is not None
# color converting
img1 = cv.cvtColor(img1, cv.COLOR_BGR2RGB)
img2 = cv.cvtColor(img2, cv.COLOR_BGR2RGB)
img3 = cv.cvtColor(img3, cv.COLOR_BGR2RGB)
img4 = cv.cvtColor(img4, cv.COLOR_BGR2RGB)
img5 = cv.cvtColor(img5, cv.COLOR_BGR2RGB)

# getting the matching pairs using SIFT algorithm
matches1to2, key1_1to2, key2_1to2 = sift_matching(img1, img2, 1, 2)
matches2to4, key1_2to4, key2_2to4 = sift_matching(img2, img4, 2, 4)
matches4to5, key1_4to5, key2_4to5 = sift_matching(img4, img5, 4, 5)

matching pairs for img1to2

matching pairs for img2to4

matching pairs for img4to5

# getting the direct matching pairs from img1 to img5
m, k1, k2 = sift_matching(img1, img5, 1, 5)

matching pairs for img1to5

# getting the homographies
H_1to2 = homoAndGraph(matches1to2, img1, img2, key1_1to2, key2_1to2)
H_2to4 = homoAndGraph(matches2to4, img2, img4, key1_2to4, key2_2to4)
H_4to5 = homoAndGraph(matches4to5, img4, img5, key1_4to5, key2_4to5)
# calculating the final homography for img1to5
H = np.linalg.inv(np.matmul(H_1to2,np.matmul(H_2to4, H_4to5)))

matches count : 1360
iteration : 566
inliers count : 1054

matches count : 1182
iteration : 566
inliers count : 554

matches count : 1678
iteration : 566
inliers count : 1170

# calculating the absolute error+
# the given homography
H_g = np.array([[6.2544644e-01, 5.7759174e-02, 2.2201217e+02],
   [2.2240536e-01,   1.1652147e+00 , -2.5605611e+01],
   [4.9212545e-04 , -3.6542424e-05   ,1.0000000e+00]])
abs_error = np.sum(abs(H_g-H)) # error

print('calculated Homography :\n', H)
print('absolute error :', abs_error)

calculated Homography :
 [[ 6.09532066e-01  4.91010556e-02  2.14775416e+02]
 [ 2.25158651e-01  1.12334404e+00 -3.11842823e+01]
 [ 4.93047329e-04 -7.17841304e-05  9.68468463e-01]]
absolute error : 12.916189792523666

# warping the img1
out_img = cv.warpPerspective(img1, np.array(H), (800, 800))
# stitching img1 onto img5
out_img[0:img5.shape[0], 0:img5.shape[1]] = img5

# implenting the results
fig, ax = plt.subplots(1, 3, figsize=(30, 10))
ax[0].imshow(img1)
ax[0].set_title('image1', color='white')
ax[0].set_xticks([]), ax[0].set_yticks([])

ax[1].imshow(img5)
ax[1].set_title('image5', color='white')
ax[1].set_xticks([]), ax[1].set_yticks([])

ax[2].imshow(out_img) 
ax[2].set_title('final output image', color='white')
ax[2].set_xticks([]), ax[2].set_yticks([])

plt.show()